Stability of cylindrical domains in phase-separating binary fluids in shear flow 
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The stability of a long cylindrical domain in a phase-separating binary fluid in an external shear 
flow is investigated by linear stability analysis. Using the coupled Cahn-Hilliard and Stokes equa- 
tions, the stability eigenvalues are derived analytically for long wavelength perturbations, for arbi- 
trary viscosity contrast between the two phases. The shear fiow is found to suppress and sometimes 
completely stabilize both the hydrodynamic Rayleigh instability and the thermodynamic instability 
of the cylinder against varicose perturbations, by mixing with nonaxisymmetric perturbations. The 
results are consistent with recent observations of a "string phase" in phase-separating fluids in shear. 

PACS numbers: 68.10.-m, 64.75.-fg, 47.20. Hw, 47.20.-k 



I. INTRODUCTION 



Phase-separating binary fluids form complex patterns 
of domains after a quench into the two-phase region of the 
phase diagram. The domain morphology is determined 
by a number of factors such as the volume fractions of the 
two phases, their viscosities, and any external forces ap- 
phed to the system ||l|,| . Of particular interest here is the 
effect of an external shear flow applied to the fluid. The 
shear flow competes with the phase-separation process, 
influencing the morphology and stability of the domains. 
Besides being a fascinating problem in nonequilibrium 
physics, this question is of practical signiflcance because 
the flnal properties of industrial materials involving bi- 
nary fluids often depend on the domain morphology. 

At late times after a quench into an unstable state, a 
phase-separating binary fluid consists of domains of the 
two phases which coarsen with time. The presence of a 
shear flow dramatically alters the kinetics of the phase 
separation. The effects of the deformation by the shear 
flow depend on the relations between the various time 
scales in the system. The characteristic time scale for the 
shear flow is just the inverse shear rate, I/7. When this 
time is shorter than the characteristic relaxation time of 
thermal concentration fluctuations tj, jr^ > 1, the sys- 
tem is in a "strong-shear" regime in which the critical 
fluctations are modifled by the flow. Conversely, in the 
"weak-shear" regime < 1 the critical fluctuations are 
not affected by the flow. A third time scale, which will 
be crucial here, is associated with the domains. Clearly 
when the growth rate of the domains, or the growth 
rate of any instabilities associated with the domains, is 
smaller than the time associated with the flow, the shear 
flow will affect the morphology and stability of the do- 
mains. 

Of interest in this paper is the competition between 
the shear flow and the coarsening process. The shear 
flow tends to deform and elongate or fragment domains, 
whereas the thermodynamics favors coarsening to larger, 
isotropic domains. This competition leads to the for- 



mation of a nonequilibrium, dynamic steady state in 
which the coarsening is stopped by the shear flow 
When the viscosities of the two phases are similar and 
one phase forms a droplet phase, the steady state con- 
sists of somewhat deformed droplets of the order of the 
Taylor breakup size R ~ o'/rjj, where a is the surface 
tension, rj the viscosity and 7 the shear rate On 
the other hand, when the two phases are both percolated 
the morphology is an anisotropic bicontinuous structure 
with apparently stable domains highly elongated along 
the flow direction. The anisotropy in these bincontinu- 
ous patterns is much larger than the aspect ratio of 2 or 
3 seen for isolated droplets Microscope observations 
have shown that the domains can be elongated into long 
cylinders ||l0|-p^. In weak shear, 7t^ < 1, these string- 
like domains still undergo frequent breakup, reconnec- 
tion, and branching, whereas in the strong-shear regime 
the system forms a "string phase" consisting of macro- 
scopically long cylindrical domains aligned with the flow 
direction. These are surprising observations, since a long 
fluid cylinder at rest is unstable against breaking up into 
spherical droplets via the Rayleigh capillary instability 
|p^ , p^ . In the situation under consideration here, the 
string is a domain of one phase immersed in the other, 
which we would also expect to be thermodynamically un- 
stable since the cylinder could lower its surface energy 
by spheroidizing. Thus the shear flow stabilizes both the 
thermodynamic instability towards phase separation and 
the hydrodynamic instability of these highly elongated 
domains. 

The goal of this paper is to explore the stabilization 
of cylindrical domains by an imposed shear flow. It is a 
sequel to my previous work on the stabilization by shear 
flow of a two-dimensioiial, lamellar domain in phase- 
separating binary fluids |15|| . In that paper it was shown 
that a lamellar domain at rest with diffuse interfaces is 
unstable towards a "varicose" instability. This instabil- 
ity is essentially a coarsening effect and depends on the 
flnite width of the interfaces; in the limit of mathemat- 
ically sharp interfaces a lamellar domain is stable. An 
external shear flow stabilizes the lamellar domain by ad- 
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vecting the top and bottom interfaces with respect to 
each other so that they no longer can maintain the ex- 
act phase relation which produces the unstable varicose 
mode. 

The instability of a quiescent cylindrical domain of one 
phase immersed in the other is somewhat different, ow- 
ing to the different dimensionality. There are two sep- 
arate forces driving instability, one hydrodynamic and 
one thermodynamic. Rayleigh jl3| and later Tomotika 
analyzed the instability of an infinitely long viscous 
fluid cylinder immersed in an immiscible fluid to axisym- 
metric varicose perturbations as shown in Fig. 0. When 
the wavelength of the varicosity is equal to or longer than 
the circumference of the cylinder the perturbation is un- 
stable and grows. This is because the higher curvature in 
the necks leads to a higher Laplace pressure there than 
in the bulges, which tends to drive fluid from the necks 
toward the bulges. Eventually the cylinder will break 
up into spherical droplets, with less total surface energy 
than the original cylinder. However, even in the absence 
of fluid motion (e.g. consider a cylindrical domain in a 
solid binary alloy) a cylindrical domain in the two-phase 
state is still unstable. This is due to the Gibbs-Thomson 
effect, in which the chemical potential depends on the 
curvature . The higher curvature at the necks will 

drive a diffusive flux towards the bulges, also leading to 
instability. Both of these mechanisms are present even in 
the limit that the interfaces are mathematically sharp. 

In this paper I perform a linear stability analysis to 
investigate the effect of an external shear flow on the sta- 
bility of a single infinitely long cylindrical domain, per- 
fectly aligned with the flow. I consider late times after a 
quench into the two-phase region, when the fluid consists 
of domains of the two phases close to their equilibrium 
concentrations, separated by well-defined interfaces. The 
formulation of the problem allows in principle for dif- 
fuse interfaces with a finite width ^ (when the viscosi- 
ties of the two phases are equal) , but in practise results 
are much easier to obtain in the limit of mathematically 
sharp interfaces and the corrections due to finite widths 
will not affect the results qualitatively. The shear rate 
7 is assumed to be small enough that the system is in 
the weak-shear regime, < 1, so that the shear fiow 
will not influence the structure of the interfaces them- 
selves. In this case the usual hydrodynamic equations for 
a phase-separating fluid are valid. As mentioned above, 
the string phase itself seems to form only in the strong- 
shear regime However, here I will be concerned not 
with the formation of the string phase but with its stabil- 
ity. The goal is to understand the mechanism by which 
the shear stabilizes these remarkably elongated domains. 
The results may also illuminate the stability of the highly 
anisotropic, bicontinuous morphologies observed in weak 
shear. I will neglect the ends of the string and also the 
possiblity that it could be inclined at a small angle to 
the flow direction. This approximation seems reason- 
able given the extraordinarily high aspect ratio observed 
for the strings and the fact that long slender drops in 



shear have a long central portion which is cylindrical and 
aligned with the flow [p^ . 

As well as shedding light on the stability of elongated 
domains in phase-separating fluids, this work encom- 
passes the problem of the effect of shear flow on the 
purely hydrodynamic viscous Rayleigh instability in im- 
miscible fluids (neglecting the thermodynamic effects). 
To my knowledge this problem has not been solved be- 
fore in the particular limit examined here. Russo and 
Steen |]l9| and Lowry and Steen [ pOpl| ] found that ax- 
ial flows can suppress capillary instabilities on cylindri- 
cal interfaces. Most other studies in the fluid dynam- 
ics literature concerning the effects of externally imposed 
flows on long fluid cylinders have been in the context of 
drop breakup. Several authors have considered the linear 
stability of an infinitely long fluid cylinder in an elon- 
gational flow p^-p5[. The flow field limits the growth 
of any disturbance to a finite value so that there is no 
true instability, and the cylinder is stabilized. However, 
some disturbances have time to grow transiently to a fi- 
nite amplitude comparable to the decreasing radius of 
the elongating cylinder, causing breakup (even though 
the disturbances do not grow exponentially). Khakhar 
and Ottino extended the analysis to general linear 
flows including shear flow, but only in the case of small 
asymmetry, when the shear part of the flow is small com- 
pared to the stretching. In this paper I will explore the 
opposite limit in which the stretching is negligible but 
the asymmetry is large. Finally, Hinch and Acrivos stud- 
ied a finite, long slender drop in shear flow |l|l. They 
find steady-state solutions for the shape of the drop at 
all shear rates, but these equilibrium solutions are un- 
stable above a critical shear rate. This is essentially due 
to the fact that the ends of the drop are not completely 
aligned with the flow, so that at sufficiently high shear 
rates the drop cannot balance the stretching of the ends 
and it extends transiently, becoming progressively thin- 
ner. This does not happen for an infinitely long cylinder, 
as considered in this paper. 

In Sec. |l| I will describe the model equations of mo- 
tion used to describe the fluid. In Sec. HI the equations 
of motion are linearized for small perturbations about a 
cylindrical domain. Approximate solutions can be found 
by writing the stability eigenvalue equation as a matrix 
equation in a truncated set of "basis" states, correspond- 
ing to different perturbation modes of the cylinder. The 
matrix elements are calculated in Sec. IV. We will see in 
Sec. ^ that the shear flow has the effect of mixing nonax- 
isymmetric disturbances with the axisymmetric varicose 
mode, leading to stabilization in some circumstances. I 
will first discuss the results for the special case in which 
the viscosities of the two phases are equal, and then gen- 
eralize to the case of arbitrary viscosity ratio. Some dis- 
cussion of the relations of th is work to experiment will 
be presented in Sections V C and 



2 



II. MODEL EQUATIONS 



I use the same equations of motion as in 1 15 1. A simple 
binary mixture can be described by one scalar order pa- 
rameter the difference in concentration between the 
two components. Since we are interested in late times 
after a temperature quench when the system consists of 
well-defined domains, the usual Ginzburg-Landau form 
for the coarse-grained free energy of a symmetrical mix- 
ture is sufficient to describe the thermodynamics: 

^Jdr (^Ik (V<f>)^ - iro*^ + iff<i>^) , (2.1) 

where rg and g are positive constants so that the fluid 
is in the two-phase region. Minimizing the homogeneous 
part of F leads to the values of the concentration in the 
two bulk phases at equilibrium: 



The fluid is assumed to be incompressible and sufficiently 
viscous that inertial effects are negligible. The equations 
of motion for the system are then the modified Cahn- 
Hilliard equation for <i>, the Stokes (creeping flow) equa- 
tion for the velocity field u, and the incomprcssibility 
condition: 



9$ 



,dF 



— = -u • V$ + MV^ — , 



dt 



, 6F 
(5$ 

= V • u. 



VP, 



(2.2) 

(2.3) 
(2.4) 



Here M is a concentration-independent mobility; v is the 
viscosity; and P is the pressure, which in general is de- 
termined by the incomprc ssib ility condition (2.4). The 
equation for the velocity (2.3) is generalized to include 
the coupling of the order parameter to the velocity field 
p6| . This term leads to a capillary force at inter face s, 
wher e gradients in $ induce fluid flow. Equations ( |2.2| )- 
(2.4) are the same as those of "model H" (without the 
thermal noise terms) used to study critical binary flu- 
ids These equations have been used extensively to 
study phase separation in binary fluids [p8| . 

Now consider a single cylindrical domain of radius TZ 
composed of say, phase a with viscosity i^*, immersed 
in an infinite region of phase (3 with viscosity i>° as il- 
lustrated in Fig. ^. The external shear flow is imposed 
along the x direction by applying a constant shear stress 
Ho far from the cylinder. Below I will allow for a finite 
width interface between the two phases only in the case 
that the viscosities are equal, i/'^ = = v] when the 
viscosities between the two phases are different I will as- 
sume the interfaces are sufficiently sharp that the viscos- 
i ty c hanges discontinuously at the interface so that Eq. 
(2.3) holds in the two different phases separately. The 



first step in a stability analysis of the cylinder is to de- 
rive the steady-state solutions to the equations of motion 
which correspond to these assumptions and to the geom- 
etry of Fig. H We therefore assume that $ is a function 
of r only and that the velocity is only nonzero in the x 
direction, u ~ u{r,9)ii, and look for ti me- independent 
solutions. The Cahn-Hilliard equation ( ^.2| ) has steady 
state solutions satisfying 



5F 



-XV $ - To* -I- .g* = M = const, 



(2.5) 



where fi is the exchange chemical potential. Using cylin- 
drical coordinates (r, 9, x) the stationary concentration 
profile 4>s(r) therefore satisfies 

(P(j)s Kd<j)s 3 , , 

1 ra<l)s + g<f>l = ^J.- (2.6) 

z r ar 



dr 



For a sufficiently large radius, (t>s will approach the profile 
for a flat interface between the two coexisting phases. 



(j)s (r) = (f>e tanh [(r - 7^)/^] , large TZ, 



(2.7) 



where the width of the interface between the two coexist- 
ing phases is the thermal correlation length ^ = y/2K/rQ 



I will assume throughout that TZ » ^, so that Eq. (2.7) 
is reasonable. Often it will be justified to further ap- 
proximate the interfacial profile by a step function so the 
interfaces are sharp, 



0,(r) «0ee[(r-7^)/^] 



(2.8) 



Note that for either interfacial profile there is a sur- 
face tension associated with the presence of the interface, 
which is just the excess free energy per unit area at the 
interface jl^ : 



K 



dr 



dr 



-K^l/^ = -ro<f>lt (2.9) 



If the viscosities of the two phases are equal, i^"^ — — 
V, applying a constant shear stress Ho far from the cylin- 
der leads to simple shear flow everywhere, with stationary 
velocity field 



Us = 7yx — 7r cos 0x, 



(2.10) 



where 7 = IIo/i' is the shear rate. More generally, for ar- 
bitrary viscosity ratio = jv" the stationary velocity 
field Us — Us {r, 6')x will have a different slope in the two 
phases. Taking the interface to be mathematically sharp 
as in Eq. (2.8), we can solve Equations ( p. 3D and ( p.4[ ) for 



Us inside and outside the cylinder separately and match 
the solutions at the interface at r = TZ. The velocity field 
must be regular at the origin and correspond to simple 
shear flow far from the cylinder, so that 

lim Us = jr cos 6x, 
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where the shear rate is defined in terms of the outer vis- 
cosity, 7 = no/i^°. Solving for Ug gives 



u^(r) 



27 



jr - 



COS 0x, r > TZ. 

(2.11) 



It is convenient to rewrite the equations in dimension- 
less form by scaling lengths by the correlation length ^, 
the concentration by its equilibrium magnitude in the 
bulk phases (j)e, and time by the characteristic diffusion 
time Tj. The velocity is scaled by the correlation length 
over the diffusion time: 



r = r 



t = t 



$ = 



u = u- 




p = p 



2Mro 



Note that the new dimensionless correlation length is 
^ = 1. In dimensionless form the equations of motion 
are now 

^ = -u- V# + [-iv^i - $ + l>^) , (2.12) 



= V^u + -V$ f -iv^$ - $ 



= V • u. 



-VP, (2.13) 



(2.14) 



The equations are characterized by a dimensionless pa- 
rameter, the rescaled viscosity 77: 



V 



Mgv _ AMroV 
K " 3ctC ' 



(2.15) 



(In the case of two different viscosites there are two di- 
mensionless parameters, jf and 77°.) In dimensionless 
form the stationary solutions corresponding to the cylin- 
drical domain in shear are 



4>s{t) ~ tanh(i — P), 



u,(f) 



25 



Sf- 



r cos ^x, 
/25P2 



SR ) — 



(2.16) 

r <R 
cos 6'x, f > R, 

(2.17) 



where the dimensionless radius of the cylinder is P = 
7^/^. The dimensionless shear rate S = -jr^ is simply the 
product of the shear rate and the diffusion time and 
thus represents a second dimensionless parameter that 
characterizes the strength of the shear flow. 



III. STABILITY ANALYSIS 

In this section I present the strategy for calculating 
the stability of the cylindrical domain. We know that 
an infinite cylinder is unstable to varicose perturbations. 
One could imagine other, nonaxisymmetric perturbations 
of the cylinder as well, such as the "undulation" mode 
shown in Fig. ||. In Sec. |^ we will find that the cylinder 
is stable to all of these perturbations. The shear flow will 
have the effect of mixing the different possible perturba- 
tions. 

Consider small perturbations about the stationary so- 
lutions found above (in the rest of the discussion I will 
drop the bars over the dimensionless variables for clar- 
ity): 



V = u u. 



(3.1) 
(3.2) 



To linear order in the small perturbations and v, the 
equations of motion are 



— = -Us{r) cos 6*— - (l>s{r)vr 



1 



(3.3) 



= V^v -f -4)'^{r) --V^ + Ws{r) 0f - - VP, (3.4) 



(3.5) 



= V-v. 



Here 77 is the appropriate viscosity for whichever region 
is under consideration, Vj. is the r component of the per- 
turbed velocity field v, and primes indicate differentia- 
tion with respect to r. Ws is a function of the stationary 
concentration profile: 



Ws{r) 



Mr) 



= -l + 3(/.2(r). 



(3.6) 



The time dependence of the pertu rbat ions is deter- 
mined by the concentration equation (B.2). The system 
is translationally invariant in the x direction, so we can 
write any perturbation as a sum over Fourier modes in 
X. Since we are interested in the growth (or damping) of 
perturbations we take 
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(j) = 0(r, 6)e 



ikx — ijjt 



v(r, 0)e 



ikx—Lut 



(3.7) 



We will be interested in long-wavelength fluctuations 
for which the dimcnsionless wave number fc << 1 (let 
K = /j/^ be the wave number and zu — lu/t^ be the 
growth (dampin g) ra te in the original variables). Substi- 
tution into Eq. (3.3) leads to an eigenvalue equation for 
the growth rate u: 



Luoir 



e)=ikus{r) cos e'0(r, 61) + (f>', {r)vr (r, 6) 



r dr 



d 



(--- 


d 




\ 2r 


dr 





-pi 002 - 



(3.8) 



A real, positive value of cj(fc) indicates stability of the 
cylinder agai nst t he perturbation. Since Vr depends on cf) 
through Eq. (3^), this eigenvalue equation is essentially 
an integro-differential equation in which the expression 
for Vr acts as an integral operator on cf). 

Eq. (3.8) cannot be solved exactly, so we need an 
approximate approach. Following the calculational ap- 
proach outlined in (ll), first consider the Cahn-Hilliard 
without the hydrodynamic terms: 



part of Eq. (3 



LU(f> = rF0, 

where we have defined the operators 
d 



ld_ 

r dr 



dr 



r2 9612 



(3.9) 



(3.10a) 





d 




1 92 


V 2r 


dr 




~ 2r2 9612 



fc2 



W,ir) 

(3.10b) 



This part of Eq. (3^) includes the dynamics of the con- 
centration field on the scale of the interface. Since S < 1 
(-yrj < 1), the shear flow acts on the scale of the do- 
mains but is not strong enough to alter the interfacial 
profile. As we will see shortly, this assumption allows us 
to find an approximate solution for (j). Let 0„ be the set 
of eigenfunctions of Eq. (3.9) and define a set of "conju- 
gate" functions by 

T<fn = (3.11) 

Then one can show that T and F are Hermitian operators 
(although their product is not) as long as the (/)„ and 0„ 
obey periodic boundary conditions or vanish at infinity. 
The eigenvalues cOn are real and the eigenfunctions and 
their conjugates are orthogonal: 

i(i>m,4'n) = / dr^^(r)(/)„(r) = for n ^ m. 



For any pair of trial functions (j)o and 0o obeying the 
same boundary conditions, there is a variational relation 
which gives an upper bound on the lowest eigenvalue uj 
US: 



< 



(3.12) 



Here the parentheses again indicate inner products. 

This variational theorem can be exploited to find solu- 
tions to Eq. ( |3.9| ) corresponding to various perturbations 
of the cylinder. Application of Eq. (3.12) requires a good 
trial function (j)Q. The smallest eigenvalues of Eq. (3.9) 
will correspond to eigenfunctions describing ^-dependent 
deformations of the cylinder, in which the interface is 
translated by a small amount but the interfacial width 
remains fixed |l^,^. Higher eigenvalues correspond to 
other deformation modes in which the structure of the 
interface changes, such as breathing modes which change 
the width of the interface. I will neglect all such modes 
here, since they are more quickly damped than the slow 
translational modes and are not important to the dynam- 
ics on the scale R. Thus, we can solv e the Cahn-Hilliard 
part of the eigenvalue equation (3.8) by using the varia- 
tional theorem with a trial function corresponding to the 
translational deformation of interest. 

The translational modes can be characterized by their 
angular dependences. Any general perturbation of the 
concentration field can be expanded as a Fourier series 
in 9: 



(3.13) 



where (pmir) is the function necessary to translate the 
interface by a small amount 5r in the f direction (the 
functio nal fo rm of the (j)m{r) will be calculated below 
in Sec. IV A). The Cahn-Hilliard part of the eigenvalue 



equation can then be rewritten as 

UJrn<p7nir) = r,„ (r)F„ (r)(/)m (r) , 



where the operators are 

_ 1 d2 Id m2 
™ ^ ^2d^ ^ 2^d^ ~^ 2^ 



1 cF 
2d^ 



1 d 
2r dr 



m 
2^ 



2^ ' 
1 



(3.14) 



(3.15a) 



+ ?>4>l-l. (3.15b) 



Each mode (/)„i(r)e™^ corresponds to a different geomet- 
rical perturbation mode of the cylinder. In the absence 
of the external shear flow, the cylindrical domain will be 
unstable to the axisymmetric, m — Q varicose mode as 
discussed in the introduction (see Fig. |^). We will see 
below that the rn — 1 mode shown in Fig. |^ is an exact 
solution to Eq. ( ^ ) for fc = 0. At fc = 0, it simply corre- 
sponds to a uniform translation of the entire cylinder and 
is thus marginally stable with eigenvalue ujiik = 0) = 0. 
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We might anticipate that the cyhnder will be stable to 
higher mod es in m as well. Note that the shear flow term 
in Eq. 



(3. 



is proportional to cos 9, so this term should 
have the effect of mixing modes with different values of 
m. 

Now consider the full eigenvalue equation, Eq. (3.8). 
We are interested in the stability of perturbations charac- 
terized by the various ^-dependent translational modes. 
The functions_fiim(r)e""^+*''^ are approximate eigenvec- 
tors of Eq. ( |3.9| ) . To solve the full equation we adopt an 
approximation similar to "tight-binding" ot k-p approxi- 
mations used in solid state physics. We assume the trans- 
lational modes are good basis states for the full problem 
and write Eq. 



as a matrix equation in this basis. 
We can truncate the matrix to only include a finite num- 
ber of states m and then diagonalize the matrix to find 
the stability eigenvalues. This is valid when the two hy- 
drodynamic terms are small enough that they only cause 
mixing among the m states included in the basis; they 
must be small relative to the distance to the next higher 
eigenvalue not included. Also, the shear flow must satisfy 
5' < 1 so that it is reasonable to only consider the trans- 
lational deformation modes. In the strong-shear regime 
S > 1 (7T^ > 1) the shear flow might couple to other 
modes that we have neglected, which alter the width or 
shape of the interfacial profile itself, since these modes 
only damp out on a time scale of ro ughl y rj. Note that 
the term containing Vr {r, 9) in Eq. ( p.^) depends on (f> 
through the hydrodynamic e qua tion (3.4). So for each 
mode m we can solve Eq. (3.4) for Vr(r,9), assuming 



that 0(r, 9) is given by the approximate basis function 
(j)m{r)e""'^ . In general the resulting velocity field can 
then also be expanded as a Fourier series. Denoting by 
v„i the solution for v obtained from substituting Eq. 
( p^ and (l){r,9) = 0„(r)e^''"*' into Eq. (Q, we can ex- 
pand 



v„i(r,6') = ^( 



(3.16) 



The 9 dependence of v will not necessarily be the same 
as that of (j)m so in general the coefficients v„m {r) in the 
sum will be nonzero even for n ^ m. 

To obtain the effective matrix equation corresponding 
to Eq. (3^) we write </)(r, 9) as a vector 



cj>{r,9) 



ei0i(r)e*« 
e202(r)e2'^ 



v 



and multiply on the left in Eq. ( p.SD by the corresponding 
conjugate vector. Here the are the amplitudes of the 
small perturbations (/>m- Recall the conjugate function </) 
is defined by = ^ so it satisfies the Poisson equation 



(3.17) 



We can expand in the same way as so that 
0(r,0,x,t)=^^„(r)e™''+''=-e-"*, 

■m 

in which case F^^m = 0m- We can easily solve for ^mif) 
using the Green's function for the Laplacian in cylindrical 
coordinates. The result is 

f-OO 

0m(r) - / r'dr'2if„(fcr>)/„(fcr<),^,„(r'), (3.18) 



where r< (r>) indicates the lesser (greater) of r and r', 
and Kr,i , Im are the modified Bessel functions. Then Eq. 
~8h becomes after multiplying on the left by 



) 



^(^^)* = (0o(r-) ^i(r)e-'« ^2(r)e 
and integrating over all 0, 



( {4>o,M^ 

{$u<t>i)Lo 






(eo\ 








£2 


I 





( (00, 'yr,OO0'J + (00, Fo0o) (00, \ikUs(i)\ + «r,O10s) 
(01, ^ikUs(t)0 +Vr,W<t''s) (01,'f^r,110s) + (01,Fl0l) 



^'r,21 







i>l, \ikUs(j)2 + Wr,120s) 
0'J (02,?^r,220;) + (02,F202) 













£2 







(3.19) 
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Solving this equation gives approximate stability eigen- 
values Lu{k) for the cylinder in the shear flow. We have 
used the definition r„i0„j = 0,„; the diagonal diffusive 
terms are then exactly the variational expression (3.12). 
The shear terms involving the stationary velocity (x S 
are completely off-diagonal, so they will indeed have the 
effect of mixing the modes. We will find in Sec. IV C that 
the off-diagonal elements involving v o nly mix modes 
that differ by ±f as written in Eq. ( 3.19| ) and that they 
are also directly proportional to the shear rate S. Thus 
in the absence of the shear flow, S = 0, the matrix is 
completely diagonal and the m modes are independent, 
with stability eigenvalues 



{4'7n,4>m) 
^m.d, S — Q. 



(3.20) 



These zero-shear stability eigenvalues are the sum of two 
terms, one due to hydrodynamic transport in the sys- 
tem and the other due to diffusive transport. Solving 
Eq. ( 3.19| ) for non-zero shear rate requires truncating the 
matrix at some mode to; since we expect only the m = 
mode to be (possibly) unstable, we might anticipate that 
only a few of the higher modes are needed to investigate 
the behavior of the to = mode under shear. 

To summarize the results of this section, the equations 
of motion were first linearized in the small perturbations 
(j) and V and expressed parametrically in terms of the 
wave number k. The perturbations of the cylinder of in- 
terest here, the translational modes, were characterized 
by their dependence on 9. A variational expression was 
introduced for the diffusive part of the problem, and the 
eigenvalue equation for the full problem was written as a 
matrix equation in the basis of the to translation modes. 
In the remaining sect ions I calculate the various matrix 
elements in Eq. (3.1£), which requires solving the hydro- 
dynamic equation for the perturbed velocity field v, and 
then solve the matrix equation itself and examine the 
results for various parameters of interest. 



IV. MATRIX ELEMENTS AND RESULTS 
WITHOUT SHEAR 



A. Diffusive contribution 



We begin by calculating the diffusive contribution to 
the matrix elements lo^a &s defined in Eq. (jp^). We 
need to determine the r dependence of the basis func- 
tions (j)rn- For Vr = Q and fc = 0, translating the en- 
tire cylindrical interface by an amount dr, 4)s{y + dr) — 
(j)s{i') + V(/)s • dr, requires adding the function 



V0S = </'sf — 0s (j') (cos 6'y + sin 6*2) 



(4.1) 



to the original interfacial profile. But the translated in- 
terface, and therefore V0s, should also be an exact solu- 
tion of (3.14). We can verify this by differentiating Eq. 
(^.6|) for the stationary solution (f)g with respect to r: 



2 dr3 



1 d2 



1 



2r dr^ 2r'^ dr 



dr 



dr 



But this is simply equivalent to 



■ m=l(i 



0, 



= 0. 
(4.2) 

(4.3) 



so (t>i{r) = 4>'s{r) is an exact eigenfunction of FiFi for 
k — 0, with eigenvalue uji — 0. We can exploit this solu- 
tion to approximate Fm4'm for general to. Since i? 3> 1 
{TZ 3> £,) and since V^cj) is only significant near the inter- 
face at r = i?, we can in general replace the term TO^/2r^ 
by TO^/2i?^ in the expression for Fm (3.15b). But then 
we would expect all to modes to have roughly the same 
radial dependence as the m = 1 mode, so we can approx- 
imate 



to 
2^ 



1 
2^ 



<f>'s- (4.4) 



For nonzero k, Fm^i't's — \k'^4''s^ so this gives 



(4.5) 



From Eq. ( 3.2C ) the diffusive part of the diagonal matrix 
element ojm is then 



1 



(4.6) 



We are assumi ng w e can approximate by the flat inter- 
face form Eq. (2.7) so 0^ = sech^ (r — i?). Since all modes 
have the same radial dependence we set (pmir) = (j)'g{r) 
in the denominator as well, which gives 



dr ( -k^r 



(0„(r),sech^(r-R)) 7o V2 
^^^)sGch4(r-R), (4.7) 

> 

r' dr' 2K„, {kr> )/„ (fcr< )sech^ (r' - R) . 

(4.8) 



The diffusive contribution to uj has thus been reduced to 
quadrature. 



Equations (4.7) and (4.8) are integrable numerically 
and are valid for diffuse interfaces of width ^. However, 
qualitatively the results are the same for sharp interfaces. 
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in which case the result can be expressed in closed form. 
The easiest way to take the sharp interface limit is to 
take 



sech^(r - R) ^2S{r- R), 



sech''(? 



R)^-S{r~R), 



(4.9) 



in all of the integrals, where 6{r — R) is the Dirac delta 
function. (This is equivalent to using the step-function 
profile Eq. ( |2.§| ).) The prefactors maintain the correct 
normalizations, so we are taking the width of the inter- 
face ^ — > while keeping the surface tension a constant. 
The delta functions then allow easy evaluation of the nec- 
essary integrals. The conjugate function 0,„ is 



ARK^{kR)I^{kr), r < R 
4RK,nikr)I„r{kR), r > R, 



(4.10) 



so the normalization integral is 



rdr 4)mir)2S{r - R) 



= 8R^KUkR)lUkR)- 



(4.11) 



Substituting into Eq. (|4.7| ) we find the diffusive term in 
the stability eigenvalue is 



rn — 1 



12R3K^{kR)lUkR)' 
In the original variables, the dispersion relation is 



(4.12) 



1 



which we note can be written independently of ^ as 
Ma K^n^ +m'^ - I 



(4.13) 



bi 7^3if™(K7^)/„,(K7^)■ 



(4.14) 



We immediately sec that the axisymmetric m — mode 
is thermodynamically unstable for kTZ < 1, i.e. for wave- 
lengths longer than the circumference of the cylinder, 
whereas the m = 1 translation mode is stable for all k 
as predicted. Fig. ^ shows the dimensionless part of the 
first three modes, 



1 



i2K„^{Kn)i„^{Kn)' 



(4.15) 



where vjm.dT^ — ^^m.d{i^TVj /Tl ■ These are the diffusive 
contributions to the stability eigenvalues in the absence 
of the shear flow. 

Fig. [s] shows th e di fference between using finite width 
interfaces in Eq. (4.7) and using the sharp interface ex- 
pression ( [4.13 ) for a cylinder of radius TZ = 4^, for the 
lowest two eigenvalues, vjq^cL and TjJi^d- The finite width 
curve is at most 15% more negative than the sharp in- 
terface curve for n7o,d, and 31% larger than the sharp 



interface curve for wi^d, at kR, = .5 (near the maximally 
unstable varicose mode). The shapes of the curves are 
qualitatively the same in both cases. As TZ increases 
and the step-function approximation for the interfacial 
profile becomes better, the difference between the curves 
grows smaller. For TZ — QS^ the differences are reduced 
to roughly 9% and 20%, respectively. Thus the quali- 
tative behavior should remain unchanged for a domain 
with sharp interfaces, but detailed comparison with ex- 
perimental or simulational data may require including 
the effect of having diffuse interfaces. 



B. Velocity field, equal viscosities 

Next consider the matrix elements involving the per- 
turbed velocity field v. This section will be limited to the 
special case in which the viscosities of the two phases are 
equal, 77* = 77° = rj. Then it will turn out that the veloc- 
ity matrix elements are completely diagonal, even with 
the shear flow, and it will be possible to obtain closed 
form expressions for Vmm(?') = '^m(r). 

Recall t hat the perturbed velocity field v satisfies 
Equations (3.4) and (3.5): 



1 



1 



= V^v - -VP + -(j)',¥(l)Y 



V 



= V-v, 



(4.16) 



(4.17) 



where F — — ^V^ -I- Wg. To solve for v, we follow a gen- 
eral procedure from Happel and Brenner | |3l| ] . Taking the 
divergence of Eq. ( 4T6) and applying the incompressibil- 
ity condition ( 4.17|) leads to a Poisson equation for the 
pressure P: 



V^P = V ■ {(j)',F(j)r). 
Expanding the pressure as P = Pm(r)e 



(4.18) 



imO+ikx ^~ut 



we 



find Pm{f) satisfies 



1 m? 1 d 

K + -PL - - k'Pra = - — (r^^F^^^), (4.19) 

the homogeneous part of which is simply the modified 
Bessel equation. Using a Green's function and requiring 
Pm(?') to be finite at r = and to vanish as r — > cx), we 
have 

Pm{r) = 

- / r'dr'K™(fcr>)/™(fcr<)-— [rV;F™0™(r')] 

r°° d /I 
dr' j-j [X™(fcr>)/„(fcr<)] i-k^r' 



2r' 



sech^(r' - R), 



(4.20) 
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where we have used Eg . (4.5). Substituting this expres- 
sion for P into Eq. ( 4.16 ) gives a vector Poisson equation 
for V. 

In cyhndrical coordinates, the r and 9 components of 
V are coupled: 



2 dve Vr IdP I 



We can solve for both together by writing 
£'m(?')e* 



Ve,m 



^im,9 ikx—uit 



'i'dr,i{r)e 



imO ikx—ut 
^ 1 



(4.21a) 
(4.21b) 



(4.22a) 
(4.22b) 



so that V0 has a phase shift relative to Vr- Eqs. ( |4.2l| ) 
then become 



n" J_ in' — (h'^ ^ n — 2™ ,9 



-P' 



(4.23a) 



9" -L 1,9' _ f Zr-2 _L ,9 2m 



rnp 



(4.23b) 



Adding these together we find 



■q rn 



rjr 



- P 



The homogeneous part of this equation is the modified 
Bessel equation, with general solutions Im+i{kr) and 
K.m+i{kr). Similarly subtracting gives 

{Qm - dm)" + i(ftn - dm)' - (k^ + ^~^) {Qm - d„,) 



— —P' I ™ p 



so now the homogeneous part of the equation has solu- 
tions Im-i{kr) and Km-i{kr). We can thus construct 
exact Green's functions for the combinations + 
and Qm — dm- The solutions to the inhomogeneous equa- 
tions are therefore 



Qm{r) +dm{r) = - / r' dr' Km+i{kr>)Im+i{kr<) -P/„(r') iPm{r') 

Jo \V vr' 

r°° /I m 1 \ 

Qm{r) - dm{r) = - / r' dr' Km-i{kr>)Im-i{kr<) -PL{r') + —Pm{r') - -(/.;F„0„(r') 
Jo \V W V J 



-(j)'^Fm4>m{r') 



(4.24a) 
(4.24b) 



All of the functions in these integrals are known, so we have reduced the solution for Vr,m{f) to quadrature. For the 
case of equal viscosities, we could therefore again integrate numerically to find the matrix elements for finite width 

interfaces. 

Instead I will proceed with the sharp interface approximation Eq. (4.9). In this limit the velocity matrix elements 
in Eq. ([3.1S|), divided by the normalization integrals, are 



^m.h 



, 4'm) Jo 



dr r(t)m{r)vr,m{r)25{r - R) = w,.,m(i?). 



(4.25) 



Thus iOm^h is give n not surprisingly by the radial velocity at the linearized position of the interface, Wr,m evaluated at 
r = R. From Eq. ( 4.20|) we find the pressure Pm is given by 



Pm{r) 



_ 2{PR^ + - 1) r kK'm{kR)Im{kr), r < R 
3R 1 kKm{kr)rm{kR), r > R. 



(4.26) 



Substituting into Eq. ( [4.24 ), solving for Vr, and performing some reductions using Bessel function identities, we find 



^m,h 



2{PR^+m^ -I) (I 



irjR 



-{m + l)Km+i{q)Im+i{q) - -^{m - l)Km-i{q)Im-i{q) 



+^q[Km{q)im+i{q) ~ Km-l{q)Im{q)] 

+ ^q'^[Km+l{q)Im+l{q) - Km-l{q)Im-l{q)][Km-l{q)Im+l{q) ~ Km+l{q)Im-l{q)] 



(4.27) 
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where we have set q = kR for convenience. These 
give the hydrodynamic part of the stabihty eigenvalues 
in the absence of the shear flow. Since I have taken the 
sharp interface Umit, this part of the calculation has been 
completely decoupl ed fro m the dynamics of the concen- 
tration field 0; Eq. (4.25) is equivalent to the usual kine- 
matic condition in hydrodynamic stability calculations. 
Thus the result Eq. ( 4.27 ) is an exact, rather than ap- 
proximate, solution for the hydrodynamic stability of a 
cylinder of one fluid immersed in a second immiscible 
one. It can therefore be compared directly with previ- 
ous results. For the varicose mode m = we obtain the 
stability eigenvalue 



2{k^R^ - 1) 



1 



Kiiq)hiq) + -qKoiq)Iiiq) 



--qK,{q)Io{q) 



(4.28) 



Putting this back in dimensional form we have 



<r{q^ - 1) 



Ki(q)h{q) + \qKo{q)h{q) 



-qK,{q)Io{q) 



(4.29) 



This expression is the same as that found by Stone 
and Brenner p2[ and may be obtained from Tomotika's 
general result for the dispersion relation of the viscous 
Rayleigh instability, in the limit of equal viscosities be- 
tween the two liquids uM. 



Note each mode in Eq. (4.27) can be written in the 
original variables as Wm^h — 'y^m,h{i'^'TV) / i^Ti., where 
VLm,h is dimensionless. This dimensionless part of the 
dispersion relation is graphed in Fig. |6| for the flrst three 
modes. Again the varicose mode is unstable for all 
K < l/TZ, whereas the higher modes are stable for all 

K. 



C. Velocity field, general viscosity ratio 

For the case of general viscosity ratio /i, it is much 
more difficult to solve the hydrodynamic equation with 
diffuse interfaces. To do so requires writing the viscosity 
as a smooth function of (/) that changes from the value 
inside the cylinder r;* to the value outside 7]° in a contin- 
uous way, which introduces extra terms into the original 
differential equation (2.3). Thus in this section I will 
start with sharp interfaces from the beginning. It is then 
more sensible to follow a different approach to solving the 
hydrodynamic equation. For sharp interfaces, the term 
coupling the total concentration $ to the total velocity 
u in Eq. ( |2.3| ) becomes 



where h is the curvature of the interface located at 
C(r) = and n is a unit vector normal to the inter- 
face. This is equivalent to the usual boundary condition 
on the jump in the normal stress across a fluid interface. 
Instead of including this coupling term in the hydrody- 
namic equation, we can just solve the usual creeping flow 
equations for the perturbed velocity field v* inside and 
v° outside the cylinder separately, and apply the appro- 
priate boundary conditions at the interface in the usual 
manner. In each region v satisfies 



V^v = -WP. = V • V, 



(4.30) 



and the pressure simply satisfies the Laplace equation. 



V^P = 0. 



(4.31) 



The location of the interface of course depends on the 
perturbation mode in question; since (t>m{f) is now a 
step function we only need know the location of the step, 
which is given by 



im6 -\-ikx 



(4.32) 



for each mode m, where e is the amplitude of the small 
perturbation. 

The solutions to Eqs. (pO|)-(|l3l|) are straightforward 
so the details are given in the appendix. The general solu- 
tions are found in terms of modified Bessel functions; the 
exact solutions are then found by applying the bound- 
ary conditions and solving numerically for the unknown 
constants. The r esults arc the coefficients v„„i(r) in the 
expansion ( 3.16| ). As shown in the appendix, the only 
nonzero coefficients are t he ones for which n — m and 
n = m ± 1 so Eq. ( 3.19| ) is tridiagonal as written. As 
was shown in the last section in Eq. ( 4.25D , the velocity 
matrix elements are then simply given by 



{<t>n,4'm.) 



a.m{r = R). 



(4.33) 



In the case of zero shear fiow, only the diagonal ma- 
trix elements are nonzero and the hydrodynamic part of 
the stability eigenvalues is therefore ujm,h = Vr,mm{R)- 
Again the results here for S" = are decoupled from 
and are therefore exact, for two immiscible liquids. 

I show in the appendix that in the absence of the shear 
fiow the velocity matrix elements and hence the eigenval- 
ues are proportional to l/ri°R, just as in the /i = 1 case. 
Figs. 0-||show the dispersion relations for the dimension- 
less parts of the three lowest modes for different values 
of fi (holding t]° constant). The results for the varicose 
m — mode are again the same as those of Tomotika 
. As the cylindrical domain becomes less viscous rel- 
ative to the background, it becomes more unstable and 
the wave number of maximum instability Kmax becomes 
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smaller. Note that the damping rate for the undulatory 
TO = 1 mode depends less strongly on ^. From Fig. jo] we 
see that the m = 2 depends only weakly on kTZ. 

Finally, in the appendix the hydrodynamic equations 
are solved analytically at fc = 0, providing a check on the 
numerical results. The stability eigenvalues at k — ar e 
diagonal and independent of the shear rate S, Eq. ( A32| ): 



i^m,h{k = 0) 



3r/°i?(/i+ 1)' 



TO > 2, 



(4.34) 



with ujQ^hik = 0) = u!i^h{k = 0) = 0. Note t hat si nce the 



only remaining off-diagonal elements in Eq. (3.19) vanish 
at fc = 0, this demonstrates that the cylinder is stable to- 
wards a;-independent deformations which would result in 
a non-circular cross-section, even in the presence of the 
shear flow. 



D. Shear flow contribution 

Finally we evaluate the off-diagonal shear flow matrix 
elements 

1 ., J0„,Ms(?-)0m) 
— tfv ~ ^ 

2 i4>n,4>m) 

where we have divided by the normalization integral. 
Since we are assuming that 4'm{i') has the same radial 
dependence for all to, we simply have in the sharp inter- 
face limit 



ik (0„,Us(?')0m) _ ik 



2 (0n,0m) 2(0„,0,„)7o 

ik 



^■U IT3\ ^^^^ 



dr rus{r)(i>n4>m. 
dr rUs{r)(f>n2S(r — R) 
(4.35) 



Thus each off-diagonal element is the same independent 

of TO. 



V. RESULTS WITH SHEAR 

With all the necessary matrix elements in hand, we 
can now solve for the stability eigenvalues in shear by 
diagonalizing Eq. ( 3.1£ ). To do so requires truncating 
the matrix at some point. Since the off-diagonal matrix 
elements are all proportional to S, we can assume that 
any perturbation which damps out more quickly than the 
time associated with the shear flow can be ignored. Thus 
we only need include modes whose damping rates are less 
than or of the order of the applied shear rate. For com- 
parison purposes, the values of the stability eigenvalues 
at fc = give good estimates of the damping rates of high 
TO modes without having to calculate the full dispersion 
relations. 



In this paper it will be sufficient to only includ e the flrst 
three modes, in which case the matrix equation ( [3.1£ ) be- 
comes a 3 X 3 secular equation for ut: 



UJQ — UJ 



ikRS 



f + iSuJol 



UJi 



ikR S 



Y + iSuJ2i 



ikB.S J 
LJ2 



iSuJi2 

- LO 



(5.1) 



where I have pulled a factor of iS out of the off-diagonal 
velocity matrix elements, Vnrn.r{R) = iSuJnm- This leads 
to a cubic characteristic equation for the stability eigen- 
values Lu{k). Since the characteristic equation has real 
coefficients, the roots can be found analytically; there 
will be either three real roots or one real root and a com- 
plex conjugate pair. 

Recall that in the absence of the shear flow, the stabil- 
ity eigenvalues are simply given by the s um o f the diffu- 
sive and hydrodynamic terms as in Eq. (^.20|): 



0J„ 



0. 



I showed in Section |IV A| and Section |IV C| that these 
two terms scale differently with the parameters of 
the system, with ujm.d = ^m,d{kR)/R^ and Wm^/i — 
^m.h{kR,fJ.)/ri°R. The relative magnitude of these two 
terms depends on the dimensionless viscosity parame- 
ter 1]° = 2£^v°/3aT^ and on the radius of the cylinder. 
For sufficiently viscous and/or small cylindrical domains, 
the diffusive term will dominate, whereas for less viscous 
and/or large cylindrical domains, the hydrodynamic term 
dominates. In the following I first examine the stability 
eigenvalues in the two extreme limits in which the hydro- 
dynamic terms or the diffusive terms can be disregarded 
entirely. I will then present some results for experimen- 
tally realistic parameter values. 



A. Stabilization of the Rayleigh instability 



First consider Eq. (5.1) in the limit that the diffusive 
terms in the diagonal elements are negligible, so that 
LOm = ^m.h- Physically this corresponds to examing 
the effect of the shear flow on the purely hydrodynamic 
Rayleigh instability. Since the diagonal matrix elements 
LUjn.h scale as l/r]°R while the off-diagonal elements do 
not depend on r]° or R, the stability eigenvalues lo will 
also scale with 1 / r fR . Denote the eigenvalues found by 
diagonalizing Eq. (|5.l|) with the shear flow present by 



^3,h{s) = ^^njj,{kR,n,s), 



(5.2) 



where the index j refers to the order of the new eigen- 
values: ujq{S) < LUi{S) < . . .. Returning to the original 
variables, this relation is simply 



(5.3) 
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Since a/v°TZ has units of inverse time, we can mea- 
sure the shear rate 7 just as weh in these units as in 
units of so for the rest of this section I will take 

S = v^n^ja. 

We start by considering the equal viscosity case, /i = 1. 
Figs. (|lo|)-(p^ show the real part of the dimensionless 
dispersion relations flj^h{kR, 1, S) for the first two stabil- 
ity eigenvalues in shear flow at various shear rates (the 
third mode was included in the calculation but is not 
a particularly interesting function of S). For some val- 
ues of K the lowest two modes uJo^h{S) and uJi^hiS) are a 
complex conjugate pair and so they show up as a single 
curve m Figs. ®-(|l2|); for the regions in k for which 
there are two curves shown, the two eigenvalues are both 
completely real. We see that at small k where there is 
a complex conjugate pair, the perturbations are travel- 
ing waves with an overall damping rate. At low shear 
rates the mode with the lowest eigenvalue LUo^h{S) is still 
unstable, but the window of wave numbers over which 
ujQ^hiS) < becomes smaller as S increases. At some crit- 
ical shear rate Sc the minimum in too^hiS) crosses zero, 
at a critical wave number Kc = kc/TZ. Above this critical 
shear rate, the instability is gone — the initially unstable 
varicose mode has been stabilized by the applied shear 
flow, by being mixed with the higher modes. 

Clearly, the critical shear rate 7c must have the same 
dependence on a, and TZ as the stability eigenvalues: 



7c 



(5.4) 



Thus the critical shear rate is a monotonically decreasing 
function of the radius of the cylinder and the magnitude 
of the outside viscosity; the smaller the cylinder or the 
less viscous the fluids are overall, the faster the growth 
of the varicose mode and so the critical shear rate must 
be faster as well for stabilization to occur. Inverting Eq. 



{ pA[) we see that at a fixed shear rate, there is a criti- 
cal radius Rc = jv° /aS{ii) above which the cylinder is 
stable and below which it is unstable. Thus if instead of 
an infinite cylinder we had a finite long cylindrical drop 
that was being stretched by the flow, initially small capil- 
lary disturbances on the drop would be suppressed by the 
flow, but as the drop thinned to a radius smaller than Rc 
the disturbances would start to grow and the drop would 
break up. Note that Sc is considerably larger than the 
magnitudes of both rio,/t(5') and Qi^h{S) at all kTZ < 1, 
so the shear rate must be faster than the rate of growth 
of the instability to stabilize it. 

The qualitative picture here is that the shear flow ad- 
vects opposite sides of the cylinder relative to each other 
so that the special axisymmetric, varicose perturbation 
no longer exists long enough to be unstable. This pic- 
ture is born out by the eigenvectors corresponding to 
the eigenvalues shown in Figs. |l0|-|l^. Fig. |l^ shows a 
cross-section of the eigenvector corresponding to the low- 
est eigenvalue at KmaxTZ = .56219 (where Kmax is the 
wave number of maximum instability in zero shear) and 
S — .1, when the lowest mode is still unstable; Fig. n4 



shows the same eigenvector at S' = .18, which is above 
the critical shear rate and therefore stable. As the shear 
rate increases, the original varicose mode becomes more 
distorted as it is mixed with the other modes. 



Next we explore the effect of the viscosity ratio be- 
tween the two phases, fj,. When the viscosities of the two 
phases are equal, above the critical shear rate the cylin- 
der is stable against perturbations at all wavelengths as 
we see in Fig. |l^. This is not the case for all /j,. For 
general /i, the shear flow does stabilize the varicose mode 
around the main instability at Kmax, but for some ^ there 
is a residual instability left at small wave numbers (long 
wavelengths). An example for fi = .25 is shown in Figs. 
[isl and ^ Fig. |l5| shows the lowest two stability eigenval- 
ues at a low shear rate, when the lowest mode ujo.hiS) is 
still unstable near the original maximally unstable wave 
number KmaxTZ — .59. We see that there is an additional 
unstable region at small k, separate from the main insta- 
bility. This long wavelength instability remains after the 
main instability has been stabilized by the shear flow, 
as shown in Fig. |l6|. The physical significance of this 
residual instability will be discussed further in Sec. VC; 



it does mean that the cylinder is still unstable to very 
long wavelength perturbations. No residual instability 
was found for .8 ^ ;< 1.0; as is either increased or de- 
creased away from this range, a small instability appears 
smoothly from k = and extends over increasingly large 
K as /i becomes correspondingly smaller than about .8 or 
larger than about 1.0. 

Nevertheless we can calculate the critical shear rate re- 
quired to stabilize the original maximally unstable mode, 
as a function of /j,. The result is shown in Fig. ^ The 
graph only includes values in the range .04 < fi < 2.4 
since for values of outside this range, the critical shear 
rate Sc becomes larger than the damping rate of the j = 2 
mode; to extend the range of /i would therefore require 
including t he j = 3 mode and higher as Sc increased. 
[From Eq. (|J|), we find n2 h{k = 0,/i = .04) = .65, 
^sAk ^0,^i^ .04) = .98, n2Ak = O,^^ = 2.4) = .20, 
and fls^k = 0,^ = 2.4) = .29 so the range in Fig. |l^ 
is reasonable.] Sc{^J■) has a minimum near fi — .5 and 
rises on either side, so that as the domain becomes either 
more or less viscous than about half the outside viscos- 
ity, it requires a higher shear rate to stabilize it. The 
rather sharp bend near ii — .1 is due to the fact that 
the maximally unstable varicose mode with growth rate 
^o,hikmaxRT l^-T S — 0) movcs to lower wave numbers as 
fi decreases (see Fig. 0), while the m = 1 mode in par- 
ticular changes less with /i (Fig. ||). The magnitudes of 
the two eigenvalues both increase as the inner viscosity 
decreases, but the growth rate of the m = mode does so 
more quickly and with a larger change in the dependence 
on kR. Thus as decreases the interaction between these 
two modes changes in such a way as to result in the fairly 
sharp increase in shear rate necessary for stabilization for 
/i < .1. 
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B. Stabilization of the thermodynamic instability 

In the last section we saw that the hydrodynamic 
Rayleigh instabihty can be partially or completely sta- 
bilized by the shear flow, depending on the viscosity con- 
trast between the two phases. The opposite limit is to 
consider what happens when the fluids are so viscous 
that the hydrodynamic terms are negligible. Then the 
diagonal elements in Eq. (5.1) are just com.d and the off- 
diagonal elements are the ones from the imposed flow Us- 
Fig. ^ and Fig. |l^ show the dimensionless parts of the 
first two modes, ^j,d{S) = ujj.d{S)R^ = mj,d{S)T^TZ^ , 
for two different shear rates at = 1. Here I have de- 
fined 5** = SR^ so that the trivial dependence on R can 
be factored out. Again as S* is increased the window of 
wave numbers for which the varicose mode is unstable 
becomes smaller. Also, the mode of maximum instabil- 
ity moves to lower k. Fig. 20 shows all three modes in 
shear for a rather high shear rate, S* = 2. The varicose 
mode has been mostly stabilized, with a small residual 
instability at small wave numbers. Including the fourth 
mode, J = 3, allows us to raise the shear rate up to some 
fraction of the magnitude of the j = A mode, Vli^.d — 64. 
(Note that including the j = 3 mode does not change the 
lower three modes e.g. in Fig. at all, since it is not 
mixed with them at low shear rates.) The instability at 
small wave numbers moves to smaller and smaller k, as S* 
is increased above 2 towards 64, but seems pinned near 
K = and never quite disappears. The changes with S* 
occur more slowly as S* is increased. 

Thus when the diffusive terms dominate the behavior, 
there is no well-defined critical shear rate for stabilization 
at all wave numbers /t. Furthermore, unlike the case of 
the Rayleigh instability analyzed above, the mode which 
is maximally unstable at 5* = does not cross the axis at 
some well-defined shear rate independently of the resid- 
ual instability at small wave number; instead the maxi- 
mally unstable mode just shifts with S* towards smaller 
K. We could however define a critical shear rate for sta- 
bilization at any given (small) wave number Kc] in this 
case the "critical" shear rate for stabilization will scale 
simply as 



oc 



1 



i?3 



(5.5) 



since only enters in the off-diagonal shear flow terms. 
Once again the critical shear rate is a decreasing function 
of the radius R, so for a given shear rate small cylinders 
will be unstable and large ones will be stable for wave 
numbers satisfying k > Kc- 



C. Relation to experiments 

In general the stability of the cylindrical domain will be 
determined by both hydrodynamic and diffusive effects, 
depending on the system parameters. In this section I 



will flrst consider the parameters relevant to Hashimoto 
et al.'s experimental system |p^ . They have studied 
phase separation under shear flow in a pseudobinary mix- 
ture of polybutadiene (PB) and polystyrene (PS) in a 
common solvent of dioctylphthalatc (DOP). They flnd a 
correlation length ^ « 0(1000) A, a surface tension on 
the order of 10~^ erg, a diffusion constant on the order of 
10-i°cm2/s, the viscosity of PB/DOP i^pb « 1.2 poise, 
and the viscosity of PS/DOB ups ~ .3 poise Q. For 
comparison with my results I will take as an example 
a viscosity ratio between the two phases of /i — .25 so 
the cylindrical domain consists of the less viscous phase. 
For the possible range of values of 77° in the experiment, 
ri° w .06 ^ .2, the hydrodynamic terms in the diago- 
nal matrix elements LUm are significantly larger than the 
diffusive terms at all reasonable values of TZ. This is 
not surprising; at the large length scale of the domains, 
7?. >> ^, we would not necessarily expect the diffusive 
terms to be important. Thus in this case the results of 
Sec. V A apply with negligible modification. The critical 



shear rate for stabilization of the cylinder at most k, for 
/i = .25 and 77° = .1, is 



1.54 



(5.6) 



Since t he t heory only applies for i? > 3 or so (for smaller 
R Eq. (2.7) is no longer a good approximation), this shear 
rate is in the weak-shear regime, < 1, and is signif- 
icantly smaller than the shear rate necessary for forma- 
tion of the string phase seen in the experiments. It is thus 
consistent that the long cylindrical domains seen exper- 
imentally are stable, since they are seen at shear rates 
that are well above the shear rate required for stabiliza- 
tion. 

Fig. |l] shows the stability eigenvalues well above the 
critical shear rate, at — 1.5. Although this shear rate 
is in the strong-shear regime > 1 for which the the- 
ory may not strictly be valid, it seems reasonable that the 
theory can be pushed into the strong shear regime at the 
large length scales 0{TZ) considered here (if in the strong- 
shear condition is replaced by the typical time scale 
for domain fluctuations, then the theory is valid here). 
For these particular parameter values, Hashimoto et al. 
found that the length of the strings seen in the experi- 
ment were on the order of 3007?,. This length could be 
explained by the resi dual instability at small wave num- 
bers discussed in Sec. V A. The wavelength of maximum 
instability in Fig. 



21 



is approximately X = 2tt/k ^ 2507?, 
so the length of the strings seen experimentally may be 
set by the residual long wavelength varicose instability in 
the shear flow. 

Next consider the case of near-critical binary fluids. 
At the critical point, /i = 1 and 77 is a universal number; 
for near-critical fluids it has the same order of magni- 
tude, so we take the critical value 77 = .7. For this value 
of 77, the diffusive terms start to become noticeable at 
small radii, although they are still signiflcantly smaller 
than the hydrodynamic terms at larger radii, e.g. i? > 6 
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or so. Fig. ^ shows the varicose m = mode without 
the shear flow, for the diffusive and hydrodynamic terms 
separately at the smallest R that is reasonable in the 
theory. The diffusive terms do change the magnitude of 
i-Uo{S) and so will have a small quantitative effect on the 
results. The stabilization by the shear flow is again very 
similar to t he p urely hydrodynamic case as illustrated in 
Figs. (|l|)-(|lf). However, even for = 1 there is now 
once again a small residual instability at small k due to 
the diffusive part of the eigenvalues that persists at high 
shear, for the small radius TZ = 3^. Note from Fig. 
that at very small k, |i37o,d| > |^^o,/i|) so it is not surpris- 
ing that the stability eigenvalues in shear resemble those 
of Fig. at small k. Finally, the critical shear rate for 
stabilization of the main instability Sc no longer scales 
exactly with at small R due to the different scaling 
of the diffusive terms {ud oc 1/R^), but the difference is 
small. Thus, for near-critical binary fluids the effects of 
diffusive transport may be observable in string-like do- 
mains for sufficiently thin strings. 

VI. DISCUSSION 

I have shown that shear flow can stabilize an isolated 
cylindrical domain in the two-phase state of a phase- 
separating binary fluid against varicose instabilities, by 
mixing the varicose mode with the other nonaxisymmet- 
ric perturbation modes of the cylinder. Essentially, the 
shear flow distorts the varicose mode by convecting one 
side of the cyhndrical interface with respect to the other, 
eliminating the special axisymmetric, pinching character 
of the varicose mode that drove the instability. Both 
the hydrodynamic Rayleigh instability and the thermo- 
dynamic, diffusive instability are suppressed by the shear 
flow, although there are residual instabilities at small 
wave numbers in the limit that the diffusive terms dom- 
inate and also for viscosity ratios /i outside the range 
•8 ^ M ^ 1.0. Other authors have considered the ef- 
fect of flows on the Rayleigh instability, but this is the 
first study focused on the effect of the nonaxisymmetric 
nature of shear flow on the Rayleigh instability. 

Comparing with the experimental results of Hashimoto 
et al, I found the mechanism presented here for stabi- 
lization of the cylindrical domains is consistent with the 
stable "string" phase seen experimentally, and that the 
lengths of the strings may be set by the residual insta- 
bilty at long wavelengths. However, it should be noted 
that the stability of a cylindrical domain in shear flow 
docs not act as a criterion for the observed relationship 
between the shear rate and the radius of the domains 
seen in the string phase. The domains in the string 
phase are formed through a dynamic process; the ob- 
served radius is not a parameter of the system (as in 



this calculation) but rather is determined through the 
self-organization process as the shear flow competes with 
coarsening in the phase-separating system. I have merely 
demonstrated a mechanism by which these macroscopi- 
cally long, cylindrical domains may be stabilized by the 
shear flow. Although a few experiments have looked at 
the breakup of the strin g-lik e domains after complete ces- 
sation of the shear flow |12| , ^5|j3^ ] , it would be interesting 
to do a careful experimental study of the shear rate at 
which the strings first begin to be unstable to see if the 
Rayleigh and/or thermodynamic varicose instabilities ex- 
plored here are the main breakup mechanisms in these 
systems. If so then the strings should be unstable below 
the critical shear rate 7c found in this paper. 

The results presented here may also shed light on why 
the shear flow can halt the phase separation and result in 
a dynamic steady state, even in the weak-shear regime, 
7rj < 1. In a concentrated phase-separating fluid when 
the two phases are both percolated so that the domains 
form a connected bicontinuous pattern, the coarsening 
is dominated by curvature effects. Qualitatively we can 
think of a piece of the interconnected structure as a cylin- 
der of fluid immersed in the other phase. This cylinder is 
susceptible to the varicose instabilities considered here, 
and particularly to the Rayleigh instability, which leads 
to breakup of the cylindrical region into spheres. Sig- 
gia |l] used this picture to explain the coarsening rates 
seen in concentrated binary fluids. Since the shear flow 
suppresses these instabilities, one might expect that it 
could stabilize an anisotropic, bicontinuous morphology 
against further coarsening. For a given shear rate, when 
the domains are relatively small they will be unstable 
and will coarsen, but once the typical length scale has 
grown to the critical radius for stabilization by the shear 
flow, TZcii), the parts of the bicontinuous structure that 
are cylindrical and aligned with the flow will no longer 
be unstable. This then provides a mechanism for the cre- 
ation of the nonequilibrium dynamic steady state seen in 
concentrated phase-separating fluids in shear flow. 
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APPENDIX: CALCULATION OF THE PERTURBED VELOCITY FIELD 



The perturbed velocity field corresponding to a perturbation of the interface given hy r — R~ ^^^"^s+ikx j^^g^ 
satisfy the hydrodynamic equations 



V^v = -VP, 



(Al) 



= V-v, 



(A2) 



V^P = 0. 



To solve for v we again follow Happel and Brenner |3l|. The velocity is expanded as in Eq. ( 3.16 ) 



v„(r,0) -^e™«v„,„(r). 



(A3) 



We will solve for each coefficient Vnm separately so we take v oc e . We start by solving Eq. (A3) for the pressure. 

Let P = P„(r)e'"*+*''^; then in general 



Pnir) 



piln{kr) r < R 
P2Kn{kr) r> R, 



(A4) 



where pi and p2 are constants. The perturbed velocity then satisfies the inhomogeneous Laplace equation ( |Al[ ), so 
the solution will consist of a general solution to the homogeneous part plus a particular solution, v = + v^. Writing 
Vx = ivx,n{'r)e™^^'^''^ , the x component of Eq. ( Al) is 



1 . , 

1'Vx,n + -^V^ 



^ ~] '^'^^'^^ = —Pn{kr) 



A particular solution is given by 



i(?Jn{kr) + ^fcr/;(fcr) r<R 
ic^K„{kr) + if-krK'Jkr) r > R, 



(A5) 



(A6) 



where tti = pi/vj^k and 1^2 = P2/'n°k. The r and components of v satisf y E q. ( 4.21 ) without the extra </> term; 
assuming they are expanded as in Eq. (4.22) the particular solutions to Eq. (Al) are 



7rifcr/4_|_i(A:r) r < R 
TT2krK'^_^_i{kr) r < R, 



(A7) 



7rifcr/4_i(fcr) r < R 
■K2krK[-^_-^^{kr) r < R 



(A8) 



(the solutions to the homogeneous part of the equation will be included in the general solution below and so are not 
needed here). The components of v'' must satisfy the continuity equation ([A^): 

dr r r ' 

This determines the constants = ni/2 and = 1^2/2. 

Next we need a general solution to the homogeneous part of Eq. ( [Al[ ). Again these are just the appropriate solutions 
of the Laplace equation: 



iciln{kr) r < R 
ic2Kn{kr) r > R, 



(A9) 
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_ J a In-i{kr) r < R 
b^Kn-iikr) r > R, 



(All) 



and once again enforcing incompressibility gives ci = (a+ + a )/2 and C2 = ^ + ^ )/2. The result is a solution 
for v„ containing six unknown constants, a^, fe^, and 7ri,2: 



if{a+ + a-)Inikr) + m [/„(fcr) + krl'^{kr)] r < R 
-{b+ + b-)Kn{kr) + Ti2 [Kn{kr) + krE'^^kr)] r > R, 



+/„+i(fcr) + a-/„_i(fcr) + TTifcr [4+1 (fcr) + /;_i(fcr)] r<R 
2 \ b+Kn+i{kr) + b-Kr,-i{kr) - vrsfcr [if;+i(fcr) + K_i(fcr)] r > i?, 



(Al2a) 



(A12b) 



r a+/„+i(fcr)-a-/„_i(fcr)+^ifcr[/;+i(fcr)-/;_i(fcr)] r < i? 
" 2 1 6+if„+i(fcr) - b-Kn-i{kr) - n^kr [K+,{kr) - K_,{kr)] r > R. 



These equations give the general solutions for the coefhcients v„,„ in Eq. ( 3.16| ) 



It now remains to apply the appropriate boundary conditions at the interface to specify the remaining constants. 
These boundary conditions will apply to the total velocity and stress fields. Letting e denote the amplitude of the 
small perturbations, recall we have 

u = u^. + ev, 
P = Ps + eP, 
n = n" + ecr, 

where Ps is the steady state pressure, 11* the steady state stress tensor, and a the perturbed stress tensor. From Eq. 



(2.11), the components of the steady state stress tensor are 



Kr^ ^^cos9, (A13) 



'— sm & r < R 



Ke={ 7,°/9S7?2 ^ (A14) 

-r]°Ssme - ^ [jr^ - SR^j sinO r> R 

(note n is symmetric). We see that IIJ,. is continuous across the interface but 11* g has a jump across the interface 
when 77* 1]°. The difference in the steady state pressure Ps across the interface is simply the Laplace pressure across 
a cylindrical interface cf/TZ, which in our dimensionless variables (the dimensionless surface tension is ct = 2/3) is 

Ps-P: = ^- (A15) 

The boundary conditions are js^ : 

• Continuity of the velocity across the interface, — u° 

• Continuity of the tangential stress across the interface, IIJ — 11° 

• Jump in the normal stress across the interface due to the mean curvature H, 11^ — 11° — 411/3 (here H is 
dimensionless) 

To apply these boundary conditions we need to evaluate the appropriate components of the stress tensor on the 
deformed interface . The location of the cylindrical interface for mode m is r = R - ee™^+*'=^ so the unit normal 
to the interface is 
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R 



ee 



imO+ikx / 



+ ikee 



imd+ikx 



(A16) 



The hydrodynamic force on the perturbed surface is F = 11 • n. To fowest order in e, the three components of F are 



Fg = ecTer ^ee 

K 



F:r - ni. 



ivad+ikxT-rs 



■ 7 

ikee 



iraO-^ikx 



P. 



The normal stress is simply 

F„ = F • n = -Ps + tdrr + 2ifcee*™''+''^^n^^ + 0{e^). 
The two tangential components of the stress can be found from Ft = F — P„n, giving 

{Ft)9 - eagr + iA:ee^'"^+*'=^n^, + 0(6^), 



Finally, the mean curvature H is 

H 



2R 



k^R^ + - 1 

i 1 



(A17) 
(A18) 

(A19) 

(A20) 

(A21) 
(A22) 

(A23) 



the first term is the steady state pressure difference given by Eq. (A15). Since the stationary velocity already 
satisfies the boundary conditions, the perturbed velocity v must satisfy them separately. Denote the difference between 
quantities inside and outside the cylinder at r = R b y [[/ ]] = — f°- Then keeping in mind Eq. (A15) and that H^^ 
is continuous across the interface, from Eqs. ( A20 )-( A22 ) the six boundary conditions become: 



[K]] = [[ve]] = [[vx]] = 0, 
, „ 2k'^R^ +m'^ -1 



im 



[Wxr]] 



R 



(A24a) 
(A24b) 
(A24c) 
(A24d) 



Since 11^^ depends on 6, for a given perturbation mode m of the cylinder the solution v„j(r, 0) given by Eq. ( 3.16| ) 
will contain more than one coefficient v„m. The perturbed stress tensor cr depends only on v so we can write 



in9-\-ikx 



(A25) 



Noting that sin 6* = (e — e )/2i and using Eq. ( [A14| ) at r = i? we can write 



/i + 1 ^ + 1 



(■m-i)e^ 



(A26) 



Subst ituting in the full sums in Eq. ( [3.16 ) and Eq. (A25) and matching terms with the same 9 dependence in Eqs. 
(A24) then gives for the boundary conditions on each coefficient v„,„, 



_ 2 PR^ + - 1 ^ 

" 3 W ™"' 

kS 



(^er,nm]] — iv' ~ Vl ~T (^n=m+l — ^n=m-l) , 

M + 1 



mS 



R{lJL + 1) 



((^n— m+1 ^n— m— l) • 



(A27a) 
(A27b) 

(A27c) 
(A27d) 
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The right-hand side of the normal st ress co nditio n Eg. ( A27b) is only nonzero for n = m, and the right-hand sides of 
the tangential stress conditions Eqs. (A27c) and ( A27d ) arc only nonzero for n = m± 1. This immediately shows that 
for a given perturbation mode m of the interface, the o nly n onzero coefficients Vnm will be those for which n — m, 
n = m + I, and n — m — 1 . The eigenvalue equation ( 3.19 ) is thus tridiagonal. To solve for e ach matrix element 
Vr,nm{r = R) (sc c Scc. 1 V C ) , wc just Substitute the appropriate general solution from Eqs. (A12) into the boundary 



conditions ( |A27| ) and solve the resulting system of algebraic equations for the unknown constants. This is best done 
numerically given the algebra involved and was solved using a standard algorithm [3^ ] . 

We can however see analytically how the matrix elements depend on 77° and R. First consider the boundary 
condition on the normal stress, Eq. (A27b). Writing out the stress tensor we have 



2r,- 



dVr 



dr 



-piIraikR) + 2r]'kv'^^^{kR) + p2K„,{kR) - 27]°kv'^°^{kR) 
3 i?2 ' 



using Eq. (A4) for the pressure. The primes indicate differentiation with respect to the argument of Q 
Dividing both sides by r]°k leaves 



(A28) 
= kR. 



2v''° 



(9) 



1 



2,-q°R 



q 



(A29) 



The left-hand side now depends only on the six integration constants and on the dimensionless parameters q and fi 
(which remain dimensionless when written in the original variables). For the tangential stress in Eq. (A27c), 



1 dvr 



dve 
dr 



V0_ 

R 



V ( ^-^Vr^nikR) + kv'g ,^{kR) - ^Vtheta,n{^R) 



kS 



1 



(f^n— m+l <^n— m — 1) ■ 



Dividing both sides by ri° and multiplying by R gives 



«"i<„(9) + IVfl^Jq) - vl^Jq) - i7rLV°^^{q) - qVg^Jq) + v^Jq) 



qS{S, 



'n— rn+l 



(A30) 



Similarly, from Eq. (A27d) we have 



Dividing by r]°k gives 



dl'x ^ dvr 
dr dx 

[b {kv'^,nikR) + ikvr^nikR))]] 
mS 



Rifx + r 



n— m+1 



(A31) 



Both Eq. ( A30| ) and Eq. (A31) depend only on the integration constants, g, /i, and 5. In calculating the d iagonal 
velocity matrix elements Vr,mm{R), all right-hand sides are zero except in the normal stress equation (A28), so we 
see that in this case the integration constants and thus the diagonal elements as well will scale as 1/rfR. For the 
off-diagonal elements the only nonzero right-hand sides are from the tangential stress conditions, so the off-diagonal 
elements will only depend on q, /j, and S, and not on R or the magnitude of the viscosities. In the absence of the 
shear flow this implies that the stability eigenvalues coming from the hydrodynamic terms scale as l/ri°R. 

Finally, the equations (A27) can be solved analytically at k = 0. In this case the ge neral solutions to the modified 
Bessel equation become r" and r"", and also the left-hand sides of Eqs. ( A27c| ) and ( A27d| ) are zero at fc = 0. This 
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simplifies the algebra considerably. Following the same procedure as outlined above, we find that the velocity matrix 
elements and therefore the stability eigenvalues at fc = are diagonal and independent of the shear rate S: 



3ri°R{fi + 1 



-, m>2, 



(A32) 



with uJo^h{k = 0) = u!iji{k = 0) = 0. 
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FIG. 5. Corrections to the eigenvalues for sharp interfaces ti7o,d (solid curve) and zzii^d (dashed curve) due to a finite width 
interface (+ symbols), for TZ = 4^. 
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FIG. 6. Stability eigenvalues flo,h (solid curve), fli^h (dashed curve), and Q2,h (dotted curve) for jj, — 1 in the absence of 
shear. 
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FIG. 7. Stability eigenvalue flo,h for S = at difFerent /i: /i = 10 (dashed curve), /i = 1 (solid curve), jj. — .1 (dash-dot 
curve), jj, = .01 (dotted curve), and /i — .001 (long dashed curve). 




FIG. 8. Stability eigenvalue fli^h for S = at different /i: fi = 10 (dashed curve), n — 1 (solid curve), jj. — .1 (dash-dot 
curve), fi = .01 (dotted curve), and /i — .001 (long dashed curve). 
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FIG. 9. Stability eigenvalue 572, h for S = at different /i: /i = 10 (dashed curve), fJ. ~ 1 (solid curve), fj. — .1 (dash-dot 
curve), fj, = .01 (dotted curve), and /i — .001 (long dashed curve). 
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FIG. 10. Stability eigenvalues _Re[f2o,h(S')] (solid line) and Re[Qi,h {S)] (dashed line) for fj, — 1 at S = .1. 
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FIG. 13. Real part ol the eigenvector corresponding to flo,h at S — .1, in the ^ = plane. 
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FIG. 14. Real part of the eigenvector corresponding to Qo,h a,t S — .18, in the 9 — plane. 
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FIG. 20. Real part of the dimensionless stability eigenvalues ^0,^(5') (solid curve), fli,d{S) (dashed curve), and fl2. 
(dotted curve) a,t S* — 2 and ^ = 1. 
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FIG. 21. Stability eigenvalues Re[vjo{S)] (solid line), Re[vji{S)] (dashed line), and Re[zu2{S)] (dotted line) for jr^ = 1.5, 
Ai = .25, rf = .1 and 7^ = 4^. 
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FIG. 22. Varicose mode: vj^^d (solid line), vj^^h (dotted line), and the total wo (dashed line) for rj = .7 and 72. = 3^ {S = 0). 
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